#include <cmath>
#include <fstream>
#include "constrained_value.hpp"
#include "message.hpp"
#include "function.hpp"

/** @brief Default constructor

    The object will be unusable if it is called with this constructor
    until it is assigned to an initialised function object.
*/
util::function::function() :
  filename(),
  strategy(none),
  samples(0),
  error(0)
{}

/** @brief Constructor
    @param filename The file name of the VLAB function file
    @param strategy The caching strategy to use
    @param samples The number of samples to use in the cache
    @param error The error term for matching x to the curve parameter

    If @a strategy is set to @c none, then the cache is not used, if
    set to @c all, then all vaules for the cache are calculated
    immediately and if set to @c lazy, the cache is filled with a
    value the first time it is used.

    If an unknown strategy is supplied, an error message is issued and
    the function is not initialised.

    If @a samples is less than two and the caching strategy is not @c
    none, an error message is issued and the function is not
    initialised.

    If @a error is not a positive number, an error message is issued
    and the function is not initialised.
*/
util::function::function(std::string filename, util::cache_strategy strategy, unsigned int samples, double error) :
  filename(filename),
  strategy(strategy),
  samples(samples),
  error(error)
{
  switch (strategy) {
  case none:
  case all:
  case lazy:
    reread();
    break;
  default:
    util::err << "util::function -- An unkown cache strategy was supplied.  The function was not initialised." << std::endl;
  }

  if (samples < 2 && strategy != none)
    util::err << "util::function -- There must be at least sample points, but it was " << samples << ".  The function was not initialised." << std::endl;

  if (error <= 0.0)
    util::err << "util::function -- The error term must be positive, but it was " << error << ".  The function was not initialised." << std::endl;
}

/** @brief Copy constructor */
util::function::function(const util::function& f) :
  filename(f.filename),
  strategy(f.strategy),
  samples(f.samples),
  error(f.error)
{
  control_points.clear();
  control_points= f.control_points;
  cache.clear();
  cache = f.cache;
}

/** @brief Evalutate the function at x
    @param x The value in the function domain to evaluate

    The function's domain is bounded to the range [0,1].  Thus @a x
    is bounded to that range.
*/
double util::function::operator()(double x) const {
  if (strategy == none)
    return eval(x);

  const double       xd   = x * double(samples);
  const unsigned int low  = (unsigned int)(xd);
  const unsigned int high = low + 1;
  const double       w    = xd - double(low);

  if (strategy == lazy) {
    if (!cache[low].first)  cache[low].second  = eval(x);
    if (!cache[high].first) cache[high].second = eval(x);
  }

  return cache[low].second * (1.0 - w) + cache[high].second * w;
}

/** @brief Reread the contents of the VLAB function file

    The file reading is designed to be as forgiving as possible, but
    if the file contents are completely non-recognisable, an error is
    issued and the funciton will not be initialised.
 */
void util::function::reread() {
  std::ifstream in(filename.c_str());

  if (!in) {
    util::err << "util::function -- The file " << filename << " could not be opened.  The function was not initialised." << std::endl;
    return;
  }

  in >> std::ws;
  if (in.eof()) {
    util::err << "util::function -- The file " << filename << " was empty.  The function was not initialised." << std::endl;
    return;
  }

  std::string buffer;
  in >> buffer >> std::ws;

  if (buffer == "range:") {
    // The file is the original file format.

    // Ignore the range, it is assumed to be [0,1].
    getline(in, buffer);
  }
  else if (buffer == "fver") {
    // The file is of a newer format

    // Check the version number.
    {
      unsigned int major = 0, minor = 0;
      in >> major >> minor >> std::ws;
      if (major != 1 || minor != 1) {
	util::err << "util::function -- The file " << filename << " is not a known VLAB function file format.  The function was not initialised." << std::endl;
	return;
      }
    }

    // Ignore the name
    in >> buffer >> buffer >> std::ws;

    // Ignore the samples
    in >> buffer >> buffer >> std::ws;

    // Ignore the flip
    in >> buffer >> buffer >> std::ws;
  }
  else {
    util::err << "util::function -- The file " << filename << " is not a known VLAB function file format.  The function was not initialised." << std::endl;
    return;
  }

  // Get the number of points
  unsigned int points = 0;
  in >> buffer >> points >> std::ws;

  for (unsigned int i = 0; i < points; ++i) {
    double x = double(), y = double();
    in >> x >> y >> std::ws;
    control_points.push_back(util::point(x, y));
  }

  reset_cache();
}

/** @brief Reconstruct the cache */
void util::function::reset_cache() {
  cache.clear();

  switch (strategy) {
  case none: break;
  case lazy:
    cache.assign(samples, std::make_pair(false, double()));
    break;
  case all:
    cache.resize(samples);
    for (unsigned int i = 0; i < samples; ++i)
      cache[i] = std::make_pair(true, eval(double(i) / double(samples)));
    break;
  }
}

/** @brief The evaluation of the function on from the control points
    @param x The value to evaluate at

    This function uses a binary search in the parameter space of the
    curve to find the curve parameter that matches to @a x.  When it
    finds a value for @x that is within error term, then it returns
    the corresponding y value of the curve.
*/
double util::function::eval(util::d01 x) const {
  double low  = 0.0;
  double high = 1.0;
  util::point test;

  do {
    double check = (low + high) * 0.5;
    test = P(check);
    if (test.x() > x)
      high = check;
    else
      low = check;
  } while (std::fabs(test.x() - x) > error);

  return test.y();
}

util::point util::function::P(double x) const {
  const int n = (unsigned int)control_points.size() - 1;
  const int t = 4;
  double u = x * (double(n - t) + 2.0);

  util::point sum;

  for (int k = 0; k <= n; k++) {
    double coeff = N(k, t, u);
    sum += control_points[k] * coeff;
  }

  return sum;
}

double util::function::N(int k, int t, double u) const {
  double res = 0.0;
  if (1==t)
    res = Nk1(k, u);
  else
    res = Nkt(k, t, u);
  return res;
}


double util::function::Nk1(int k, double u) const {
  if (Uk(k) <= u)
    {
      if (u < Uk(k + 1))
	return 1.0;
    }
  return 0.0;
}

double util::function::Nkt(int k, int t, double u) const {
  double sum = 0.0;
  int div = Uk(k + t - 1) - Uk(k);
  if (div)
    sum = (u - Uk(k)) / div * N(k, t - 1, u);

  div = Uk(k + t) - Uk(k + 1);
  if (0 != div)
    sum += (Uk(k + t) - u) / div * N(k + 1, t - 1, u);

  return sum;
}

int util::function::Uk(int j) const {
  const int n = (unsigned int)control_points.size() - 1;
  const int t = 4;
  if (j < t)
    return 0;
  if (j > n)
    return n - t + 2;
  return j - t + 1;
}
